Predykcja opadów atmosferycznych¶
W tej ramce danych zajmiemy się próbą predykcji opadów atmosferycznych dla dnia następnego. Operować będę na ramce danych udostępnionej pod następującym adresem: https://www.kaggle.com/jsphyg/weather-dataset-rattle-package?select=weatherAUS.csv
Zacznijmy klasycznie od wgrania wszystkich potrzebnych bibliotek oraz załadowania ramki danych.
import numpy as np
import pandas as pd
import keras
import seaborn as sns
from matplotlib import pyplot as plt
df_AU = pd.read_csv("australia.csv")
#Ostrzeżenia dotyczące pakietów, nie przejmujmy się nimi :D
import warnings
warnings.filterwarnings("ignore")
Zobaczmy jak prezentuje się nasza ramka danych, pierw w wersji "surowej":
df_AU.head(5)
#Sprawdźmy czy nie zostaliśmy okłamani z brakiem danych
df_AU.info()
Wspaniałe więści! Wszystkie dane są na miejscu, co więcej wszysktie dane są danymi ciągłymi, co dość ułatwi nam analizę. Jak widzimy nasza ramka posiada 56420 wierszy, jest to informacja newralgiczna przy wyborze zbiorów treningowego i testowego.
df_AU.shape
Zanim przejdę do operacji na naszej ramce danych, choć zostało to już dokonane przez prowadzącego w poleceniu zadania, prztocze ze względów dostępności do informacji, informacje a propos każdej z kolumn naszego zbioru danych:
MinTemp - Minimalna temperatura [C]
MaxTemp - Maksymalna temperatura [C]
Rainfall - Suma opadów [mm]
Evaporation - Miara odparowywania [mm]
Sunshine - Suma czasu nasłonecznienia [h]
WindGustSpeed - Najwyższa prędkość wiatru [km/h]
WindSpeed9am - Prędkość wiatru o 9:00 [km/h]
WindSpeed3pm - Prędkość wiatru o 15:00 [km/h]
Humidity9am - Wilgotność o 9:00 [%]
Humidity3pm - Wilgotność o 15:00 [%]
Pressure9am - Ciśnienie atmosferyczne o 9:00 [hPa]
Pressure3pm - Ciśnienie atmosferyczne o 15:00 [hPa]
Cloud9am - Zachmurzenie o 9:00 [skala: 0 - słońce, 8 - całkowite zachmurzenie]
Cloud3pm - Zachmurzenie o 15:00 [skala: 0 - słońce, 8 - całkowite zachmurzenie]
Temp9am - Temperatura o 9:00 [C]
Temp3pm - Temperatura o 15:00 [C]
RainToday - Czy dzisiaj padał deszcz [0 - nie, 1 - tak]
Zmienna celu: RainTomorrow - Czy jutro będzie padał deszcz [0 - nie, 1 - tak]
Zacznijmy pierw od streszczonej eksploracji danych by zapoznać się potencjalnymi zależnościami.
df_AU.hist(bins = 80,figsize = (20,18))
Widzimy, że rozkłady są w większoci przypadków rozkładami Gaussa, dla zmiennych ciągłych, jednakże obserwujemy także pewne nietypowe "anomalie" w niektórych rozkładach, jak na przykład dla rozkładu zmiennej wilgotności, gdzie niekóre z wartości są bardzo nietpowo liczniejsze niż inne.
Ok, przyjrzyjmy się pairplotom i na ich podstawie wywnioskujmy, które zmienne są ze sobą skorelowane.
sns.pairplot(df_AU)
Tak, wiem, jest to całkowicie nieczytelne dla jednostki ludzkiej, zatem wybierzmy te wykresy, które demonstrują w czysty sposób pewną zależność.
sns.pairplot(df_AU,y_vars = 'MinTemp', x_vars = df_AU[['MaxTemp','Temp3pm','Temp9am']])
sns.pairplot(df_AU,y_vars = 'MaxTemp', x_vars = df_AU[['MinTemp','Temp3pm','Temp9am']])
sns.pairplot(df_AU,y_vars = 'MaxTemp', x_vars = df_AU[['MinTemp','Temp3pm','Temp9am']])
sns.pairplot(df_AU,y_vars = 'Pressure9am', x_vars = df_AU[['Pressure3pm']])
Na końcu utworzmy mapę korelacji by streścić nasze zależności:
Korelacja = df_AU.corr()
sns.heatmap(Korelacja,annot = False)
plt.show()
"Przycinanie" ramki, wstępne podejście do modelowania¶
Zanim przejdziemy do wyboru modeli uczenia maszynowego oraz ich nauki, zacznijmy od obróbki naszej ramki danych w celu ułatwienia modelom operacji na niej.
Jak już wiemy, nasza ramka nie zawiera wartości kategorycznych, a zatem w obróbce ramki zaoszczędzi nam to na treści. To samo tyczy się wartości 'null' i 'NaN', które w tabeli nie występują.
Wybrane modele wykorzystam na dwóch ramkach, na "czystej" oraz na obrobionej, będzie to zarazem test inżynieri cech, sprawdzający czy faktycznie operacje na ramce przyspieszyły, bądź nawet zmieniły (w zamyśle na lepsze) wyniki modelów.
Outliery¶
Zacznijmy od wyleminowania, niektórych outlierów. Wierszów z outlierami pozbędziemy się jedynie z "najbardziej skrajnych" kolumn, iż, jak wiemy, wartości skrajne zaburzają pracę modelu machine learningowego, jednakże dla tabelki 18-kolumnowej, w najgorszym przypdaku, moglibyśmy się pozbyć nawet 80% danych początkowych, co z pewnością nie usprawniłoby żadnego modelu.
Patrząc na histogramy naszych cech, w oczy rzuca się kilka kolumn. Mianowicie "Sunshine", "Evaporation" oraz "Rainfall", jednakże ze względu na występowanie dość znaczącej ilości wartości odstających, w szczególności outlierów minimalnych, może to przechowywać pewne infromacje dotyczące zachodzącego zjawiska, dlatego zdecydowałem się o usunięcie jedynie maksymalnych outlierów dla tych dwóch kolumn.
Rzućmy okiem na boxploty reprezentujące wyżej omawiane zmienne.
plt.figure(figsize=(12, 4), dpi=80)
plt.subplot(131)
sns.boxplot(df_AU['Sunshine'])
plt.subplot(132)
sns.boxplot(df_AU['Evaporation'])
plt.subplot(133)
sns.boxplot(df_AU['Rainfall'])
plt.show()
Jak widzimy, Evaporation oraz Rainfall są wręcz praktycznie outlierami maksymalnymi.
Na poniższych wynikach, widać, że usunięcie kolumny "Rainfall" z naszej analizy może być uzasadnione, iż około 65% wartości okupuje wartość 0, a 85% pierwsze cztery wartości, które nie różnią się od siebie znacząco.
print(df_AU['Rainfall'].value_counts()[0])
print(df_AU['Rainfall'].value_counts()[0:3].sum())
print(df_AU['Rainfall'].value_counts().sum())
#Stwórzmy kopie naszej "ramki" i przytnijmy ją
df_prep = df_AU.copy()
df_prep = df_prep.drop(['Rainfall'],axis = 1)
upper_lim = df_prep['Evaporation'].quantile(.90)
#Stosuje mniejsze przycinanie, ze względu na mniejszą ilość outlierów w tym przypadku
upper_lim_1 = df_prep['Sunshine'].quantile(.95)
df_prep = df_prep[(df_prep['Evaporation'] < upper_lim)]
df_prep = df_prep[(df_prep['Sunshine'] < upper_lim_1)]
Z outlierów maksymalnych przytnę jeszcze 5% maksymalnych wartości z kolumn "WindSpeed9am" oraz "WindSpeed3pm". Wybrałem, outliery maksymalne ze względu na "interpretację" oraz rozkład danych zmiennych. Jak wiemy stan braku wiatru jest o wiele częstszy niż występowanie dość szybkiego oraz porwistego wiatru, zatem z tych dwóch kolumn pozbędę się, jak już wspomniałem, tylko po 5% outlierów maksymalnych.
upper_lim_2 = df_prep['WindSpeed9am'].quantile(.95)
upper_lim_3 = df_prep['WindSpeed3pm'].quantile(.95)
df_prep = df_prep[(df_prep['WindSpeed9am'] < upper_lim_2)]
df_prep = df_prep[(df_prep['WindSpeed3pm'] < upper_lim_3)]
#Ile danych nam pozostało?
df_prep.shape
Skalowanie¶
Zastosujmy dla naszego modelu skalowanie. Jak wiemy niektóre algorytmy uczenia maszynowego radzą sobie "niedobrze" z danymi o różnej skali wielkości, tudzież mogą nastąpić problemy ze zbieżnością w np. gradient descencie. Jak wiemy istnieje kilka sposobów radzenia sobie z powyższym problemem. Jak wiemy istnieją dwa główne idee skalowania danych: standaryzajca i normalizacja. W tej analizie zastosuje normalizacje, gdyż sprawdza się lepiej dla różnych rozkładów, podczas gdy standaryzacja jest najbardziej optymalna dla rozkładów Gaussa.
from sklearn.preprocessing import MinMaxScaler
scaling = MinMaxScaler()
df_prep[:] = scaling.fit_transform(df_prep[:])
#Od razu lepiej
df_prep
Przygotowanie zbiorów treningowych oraz testowy¶
Jak wedle w poleceniu podzielę teraz obydwa ("surowy" oraz preprocesowany zbiór danych) na osobne zbiory treningowe oraz testowe.
from sklearn.model_selection import train_test_split
X_train_raw, X_test_raw, y_train_raw, y_test_raw = train_test_split(df_AU.iloc[:, :-1], df_AU.iloc[:,-1], random_state=24, test_size=0.4)
X_train_prep, X_test_prep, y_train_prep, y_test_prep = train_test_split(df_prep.iloc[:, :-1], df_prep.iloc[:,-1], random_state=24, test_size=0.4)
Wybór klasyfikatorów¶
Nauczmy trzy wybrane klasyfikatory.
Drzewo decyzyjne¶
Zacznijmy od wzięcia pod uwagę drzewa decyzyjnego, jako iż jest dość prostym modelem machine-learningowym.
from sklearn.tree import DecisionTreeClassifier
Tree_raw = DecisionTreeClassifier(criterion='gini', splitter='best',random_state=1337,max_depth = 5)
Tree_prep = DecisionTreeClassifier(criterion='gini', splitter='best',random_state=1337,max_depth = 5)
Tree_raw.fit(X_train_raw,y_train_raw)
Tree_prep.fit(X_train_prep,y_train_prep)
Wykorzystam funkcję, które wynonałem wraz z partnerem w projekcie na milestone II, gdyż będą niezwykle przydatne. Musiałem uwzględnić max_depth=5, gdyż w przeciwnym przypadku powstawało 31-poziomowe drzewo, przeuczające się, jednakże w analizie skupie się na parametrze "criterion".
from sklearn.metrics import plot_confusion_matrix, plot_roc_curve, plot_precision_recall_curve
def print_metrics_prep(model):
print(f"Accuracy score (train): {model.score(X_train_prep, y_train_prep)}")
print(f"Accuracy score (test): {model.score(X_test_prep, y_test_prep)}")
print(f"Tree depth: {model.get_depth()}")
print(f"Leaf count: {model.get_n_leaves()}")
def plot_scores_prep(model):
fig, axs = plt.subplots(2, 2)
plot_confusion_matrix(
model, X_train_prep, y_train_prep,
values_format="d",
ax=axs[0, 0]
)
axs[0, 0].set_title("Train data")
plot_confusion_matrix(
model, X_test_prep, y_test_prep,
values_format="d",
ax=axs[0, 1]
)
axs[0, 1].set_title("Test data")
plot_roc_curve(model, X_test_prep, y_test_prep, ax=axs[1, 0])
axs[1,0].set_title("ROC curve")
plot_precision_recall_curve(model, X_test_prep, y_test_prep, ax=axs[1, 1])
axs[1, 1].set_title("Precision recall curve")
axs[1, 1].set_ylim(0, 1)
fig.set_size_inches(15, 15)
fig.show()
def print_metrics_raw(model):
print(f"Accuracy score (train): {model.score(X_train_raw, y_train_raw)}")
print(f"Accuracy score (test): {model.score(X_test_raw, y_test_raw)}")
print(f"Tree depth: {model.get_depth()}")
print(f"Leaf count: {model.get_n_leaves()}")
def plot_scores_raw(model):
fig, axs = plt.subplots(2, 2)
plot_confusion_matrix(
model, X_train_raw, y_train_raw,
values_format="d",
ax=axs[0, 0]
)
axs[0, 0].set_title("Train data")
plot_confusion_matrix(
model, X_test_raw, y_test_raw,
values_format="d",
ax=axs[0, 1]
)
axs[0, 1].set_title("Test data")
plot_roc_curve(model, X_test_raw, y_test_raw, ax=axs[1, 0])
axs[1,0].set_title("ROC curve")
plot_precision_recall_curve(model, X_test_raw, y_test_raw, ax=axs[1, 1])
axs[1, 1].set_title("Precision recall curve")
axs[1, 1].set_ylim(0, 1)
fig.set_size_inches(15, 15)
fig.show()
print_metrics_raw(Tree_raw)
plot_scores_raw(Tree_raw)
#Praktycznie ta sama Accuracy :(
print_metrics_prep(Tree_prep)
plot_scores_prep(Tree_prep)
Ok teraz zmieńmy "criterion" z 'gini' na 'entropy'.
Tree_raw = DecisionTreeClassifier(criterion='entropy', splitter='best',random_state=1337,max_depth = 5)
Tree_prep = DecisionTreeClassifier(criterion='entropy', splitter='best',random_state=1337,max_depth = 5)
Tree_raw.fit(X_train_raw,y_train_raw)
Tree_prep.fit(X_train_prep,y_train_prep)
print_metrics_raw(Tree_raw)
plot_scores_raw(Tree_raw)
print_metrics_prep(Tree_prep)
plot_scores_prep(Tree_prep)
Widzimy, iż entropia radzi sobie odrobine gorzej niż gini, jednakże accuracy, ROC curve oraz precission recall curve są praktycznie identyczne dla odbydwu przypadków.
Las losowy¶
Skoro zaczeliśmy od drzewa, to teraz przejdźmy do lasu, apriori zakładam, iż na "logikę" las powinien zwrócić lepsze wyniki niż pojedyńcze drzewo. Analizowanym parametrem będzie "n_estimators". Ponownie musiałem ustalić maksymalną głębokość, tym razem większą, gdyż powstawało zbyt dużo przeuczonych drzew w lesie.
from sklearn.ensemble import RandomForestClassifier
forest_raw = RandomForestClassifier(n_estimators= 10,max_depth=7)
forest_prep = RandomForestClassifier(n_estimators= 10,max_depth=7)
forest_raw.fit(X_train_raw,y_train_raw)
forest_prep.fit(X_train_prep,y_train_prep)
print(f"Accuracy score (train): {forest_raw.score(X_train_raw, y_train_raw)}")
print(f"Accuracy score (test): {forest_raw.score(X_test_raw, y_test_raw)}")
plot_scores_raw(forest_raw)
#Ponownie niepokojąca obserwacja zwiększenia Accuracy w zbiorze treningowym i zmniejszenia w testowym
print(f"Accuracy score (train): {forest_prep.score(X_train_prep, y_train_prep)}")
print(f"Accuracy score (test): {forest_prep.score(X_test_prep, y_test_prep)}")
plot_scores_prep(forest_prep)
Dobrze, teraz zmieńmy ilość estymatorów z 10 na 50, większa liczba drzew powinna mieć statystycznie częściej racje, lecz zobaczymy.
forest_raw = RandomForestClassifier(n_estimators= 50,max_depth=7)
forest_prep = RandomForestClassifier(n_estimators= 50,max_depth=7)
forest_raw.fit(X_train_raw,y_train_raw)
forest_prep.fit(X_train_prep,y_train_prep)
print(f"Accuracy score (train): {forest_raw.score(X_train_raw, y_train_raw)}")
print(f"Accuracy score (test): {forest_raw.score(X_test_raw, y_test_raw)}")
plot_scores_raw(forest_raw)
print(f"Accuracy score (train): {forest_prep.score(X_train_prep, y_train_prep)}")
print(f"Accuracy score (test): {forest_prep.score(X_test_prep, y_test_prep)}")
plot_scores_prep(forest_prep)
Jak widzimy, 50 drzew nie wybiera wyników dużo lepiej niż 10, z zatem nie zaprzestańmy na zwiększaniu ilości drzew w lesie i sprawdźmy jeszcze jak poradzą sobie trzy drzewa w komitecie.
forest_raw = RandomForestClassifier(n_estimators= 3,max_depth=7)
forest_prep = RandomForestClassifier(n_estimators= 3,max_depth=7)
forest_raw.fit(X_train_raw,y_train_raw)
forest_prep.fit(X_train_prep,y_train_prep)
print(f"Accuracy score (train): {forest_raw.score(X_train_raw, y_train_raw)}")
print(f"Accuracy score (test): {forest_raw.score(X_test_raw, y_test_raw)}")
plot_scores_raw(forest_raw)
print(f"Accuracy score (train): {forest_prep.score(X_train_prep, y_train_prep)}")
print(f"Accuracy score (test): {forest_prep.score(X_test_prep, y_test_prep)}")
plot_scores_prep(forest_prep)
Konkluzje? Pojedyńcze drzewo wybiera wynik dostatecznie precyzyjnie by na własną ręke stwierdzić wynik z dużym prawdopodobieństwem, a zatem większa ilość drzew w komitecie nie poprawi tego wyniku znacząco.
Logistic regression¶
Jako ostatni klasyfikator wykorzystam regresję logistyczną, coś niezwiązanego z drzewami, jako parametr wybrałem "solver", gdyż jest on najbardziej newralgiczny w podejmowaniu decyzji.
from sklearn.linear_model import LogisticRegression
#Domyślny parametr lbfgs
regression_raw = LogisticRegression(random_state=1, solver ='lbfgs')
regression_prep = LogisticRegression(random_state=1, solver ='lbfgs')
regression_raw.fit(X_train_raw,y_train_raw)
regression_prep.fit(X_train_prep,y_train_prep)
print(f"Accuracy score (train): {regression_raw.score(X_train_raw, y_train_raw)}")
print(f"Accuracy score (test): {regression_raw.score(X_test_raw, y_test_raw)}")
plot_scores_raw(regression_raw)
print(f"Accuracy score (train): {regression_prep.score(X_train_prep, y_train_prep)}")
print(f"Accuracy score (test): {regression_prep.score(X_test_prep, y_test_prep)}")
plot_scores_prep(regression_prep)
#Liblinear jest dobry dla niewielkich ramek danych, więc powinien sprawować się gorzej
regression_raw = LogisticRegression(random_state=1, solver ='liblinear')
regression_prep = LogisticRegression(random_state=1, solver ='liblinear')
regression_raw.fit(X_train_raw,y_train_raw)
regression_prep.fit(X_train_prep,y_train_prep)
print(f"Accuracy score (train): {regression_raw.score(X_train_raw, y_train_raw)}")
print(f"Accuracy score (test): {regression_raw.score(X_test_raw, y_test_raw)}")
plot_scores_raw(regression_raw)
print(f"Accuracy score (train): {regression_prep.score(X_train_prep, y_train_prep)}")
print(f"Accuracy score (test): {regression_prep.score(X_test_prep, y_test_prep)}")
plot_scores_prep(regression_prep)
Jednak różnice są niewielkie :(
#Sag oraz saga są podobnymi parametrami, stosowanymi dla dużych ramek, więc dla skrócenia raportu, użyje wyłącznie sag
regression_raw = LogisticRegression(random_state=1, solver ='sag')
regression_prep = LogisticRegression(random_state=1, solver ='sag')
regression_raw.fit(X_train_raw,y_train_raw)
regression_prep.fit(X_train_prep,y_train_prep)
print(f"Accuracy score (train): {regression_raw.score(X_train_raw, y_train_raw)}")
print(f"Accuracy score (test): {regression_raw.score(X_test_raw, y_test_raw)}")
plot_scores_raw(regression_raw)
print(f"Accuracy score (train): {regression_prep.score(X_train_prep, y_train_prep)}")
print(f"Accuracy score (test): {regression_prep.score(X_test_prep, y_test_prep)}")
plot_scores_prep(regression_prep)
#Pozostał nam wyłącznie do sprawdzenia newton-cg, może tym razem różnica będzie zauważalna?
regression_raw = LogisticRegression(random_state=1, solver ='newton-cg')
regression_prep = LogisticRegression(random_state=1, solver ='newton-cg')
regression_raw.fit(X_train_raw,y_train_raw)
regression_prep.fit(X_train_prep,y_train_prep)
print(f"Accuracy score (train): {regression_raw.score(X_train_raw, y_train_raw)}")
print(f"Accuracy score (test): {regression_raw.score(X_test_raw, y_test_raw)}")
plot_scores_raw(regression_raw)
print(f"Accuracy score (train): {regression_prep.score(X_train_prep, y_train_prep)}")
print(f"Accuracy score (test): {regression_prep.score(X_test_prep, y_test_prep)}")
plot_scores_prep(regression_prep)
Finalne konkluzje¶
Z przykrością muszę stwierdzić, iż preprocessing nie zmienił dużo (i to niestety na gorsze), więc czasami, jak zobaczyłem, większa ilość danych bez usuwania outlierów, daje lepsze wyniki. W moim przypadku zostało usunięte około 20% danych, jak widać za dużo, przy w miare średniej wielkości ramce danych z sporawą ilością parametrów (kolumn).
W kwestii funkcjonalności klasyfikatorów, okazały się być niezwykle ze sobą zgodne i średni wynik accuracy wyniósł 85%, czyli dość porządnie.
W kwestii wyboru miary oceny, zawszę (w projektach, w pracach domowych) stosowałem accuracy jako wiodący wyznacznik gdyż pozakuje nam jak bardzo "świadomie" model stwierdził podany wynik, co jest, dla mnie, najlepszą miarą wskazania jego poprawności, chociaż mogę być odrobinę stronniczy.